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Summary 

Simple, mixed finite element models are devel- 
oped for the free-vibration analysis of curved, thin- 
walled beams with arbitrary open cross section. The 
analytical formulation is based on a linear, Vlasov- 
type, thin- walled beam theory with the effects of 
flexural-torsional coupling, transverse shear deforma- 
tion, and rotary inertia included. The fundamen- 
tal unknowns consist of seven internal forces and 
seven generalized displacements of the beam. The 
element characteristic arrays are obtained by using a 
modified form of the Hellinger-Reissner mixed vari- 
ational principle. The modification consists of aug- 
menting the functional of that principle by two terms: 
(1) the Lagrange multiplier associated with the con- 
stant condition relating the rotation of the cross sec- 
tion and the twist degrees of freedom, and (2) a reg- 
ularization term that is quadratic in the Lagrange 
multiplier. Only C° continuity is required for the 
generalized displacements. The internal forces and 
the Lagrange multiplier are allowed to be discontin- 
uous at interelement boundaries. 

Numerical results are presented to demonstrate 
the high accuracy and effectiveness of the elements 
developed. The standard of comparison is taken to 
be the solutions obtained by using two-dimensional 
plate/shell models for the beams. 

Introduction 

Since the development of the comprehensive the- 
ory of combined torsion and bending of thin-walled 
beams by Vlasov in the 1930 5 s (see Vlasov 1961), 
extensive literature has been devoted to the appli- 
cation and adaptation of the theory to a variety 
of equilibrium, stability, and vibration problems of 
thin- walled structures. Reviews of the many con- 
tributions on this subject are given in a number of 
survey papers (see, for example, Nowinski 1966 and 
Panovko and Beilin 1969) and monographs (Gjelsvik 
1981; Kollbrunner and Basler 1969; Vlasov 1961; 
and Zbirohowski-Koscia 1967). Several finite element 
models, based on Vlasov’s theory, have been pro- 
posed in the literature for the equilibrium, stability, 
and vibration analysis of thin-walled beams. (See At- 
tard 1986; Barsoum and Gallagher 1970; Epstein and 
Murray 1976; Friberg 1985; Kiss 1986; Krajcinovic 
1969; Mei 1970; Rozmarynowski and Szymczak 1984; 
Tralli 1986; and Wekezer 1987.) Also, a number of 
commercial programs have thin-walled beam mod- 
els in their element libraries. (See, for example, 
Kiss 1986.) However, except for the hybrid model 
and the modified displacement model described in 
Kiss (1986) and Tralli (1986), all the other finite 
element models developed are based on a standard 


displacement formulation that requires C 1 continu- 
ity for the torsional degrees of freedom. Despite the 
documented advantages of mixed formulations in fi- 
nite element modeling, to date no publications exist 
in which mixed finite elements are used for modeling 
thin-walled beams. The present study attempts to 
fill this void. Specifically, the objectives of this pa- 
per are: (1) to present simple and efficient mixed fi- 
nite element models for the free-vibration analysis of 
curved thin-walled beams with arbitrary open cross 
section, and (2) to demonstrate the effectiveness of 
these elements by means of a numerical comparison 
with results of two-dimensional models. 

To sharpen the focus of the study, only prismatic 
thin- walled beams with open cross section are consid- 
ered. The wall thickness is assumed to be constant, 
and the material is linearly elastic and isotropic. The 
analysis can be extended to anisotropic beams with 
variable geometric and material characteristics. 


Symbols 


A 

cross-sectional area 

Ay,Ag 

effective shear areas in y- and 
^-directions, respectively 

B u 

bimoment 

b 

flange (or web) dimension 

E 

Young’s modulus 

[F] 

linear flexibility matrix for an 
individual element 

G 

shear modulus 

{H} 

vector of internal force 
parameters 

Iy,Iz,Iyz 

second moments of cross 
section (moments and product 
of inertia) 

Itjj, lu iy> Iu>z 

second sectorial moments of 
cross section 

J 

Saint- Venant torsion constant 
of cross section 

K 

kinetic energy 

L 

beam length 

l 

length of individual finite 
element 

m 

consistent mass matrix for 
individual element 


bending and twisting moments 

m 

number of displacement nodes 
in element 


N 

shape functions used for 
approximating generalized 
displacements 

N 

shape functions used for 
approximating internal forces 
and Lagrange multiplier 

N x 

axial force 

[P},[Q] 

elemental matrices associated 
with constraint condition and 
regularization term in the 
functional, respectively 

Qy > Qz 

transverse shear forces 

R 

radius of curvature of center- 
line of beam 

[5] 

strain-displacement matrix for 
individual element 


first sectorial moment of cross 
section 

s 

number of parameters used in 
approximating each of internal 
forces and Lagrange multiplier 
within individual elements 

t 

wall thickness 

U 

total strain energy 

u c 

complementary strain energy 

U n 

strain energy associated with 
normal stresses (N x , M y , M z , 
and B w ) 

U u 

strain energy associated with 
torsional shear stresses ( M t ) 

u,v,w 

displacement components in 
coordinate directions 

O 

$ 

O 

o 

9 

axial and transverse displace- 
ments at y = z = 0 

W 

vector of nodal displacements 

x,y,z 

orthogonal coordinate system 
with x passing through cen- 
troids of cross sections 

7x y ? 7x 2 5 

/-yO 

fxy> ixz 

transverse shear strains (see 
eqs. (A2), appendix A) 

e 

penalty parameter 


extensional strain 

c x 

extensional strain of centerline 
of beam 


c 

contour coordinate used 
in evaluating sectorial 
coordinate 

9 ° 

rate of twist of beam 

*.0 *.0 *.0 
"'y ? "'z * "'t 

curvature changes and twist of 
beam 

A 

Lagrange multiplier 

{A} 

vector of Lagrange multiplier 
parameters 

P 

mass density of material 

°x 5 CTxy 5 ®xz 

normal and shearing stresses 
on cross section 

V 

Poisson’s ratio of the material 


rotation components 


strain parameter 

fi 

frequency of vibration, Hz 

w 

sectorial coordinate (warping 
of cross section for a unit rate 
of twist) 

d 

= d/dx 

Range of indices: 

hJ 

1 to s 

if j 

1 to m 

Subscripts: 


2d 

two-dimensional model result 

s 

shear center 

Superscript: 


t 

matrix transposition 


Mathematical Formulation 

The analytical formulation is based on a linear 
form of Vlasov’s thin- walled beam theory with the ef- 
fects of flexural-torsional coupling, transverse shear 
deformation, and rotary inertia included. A mixed 
formulation is used in which the fundamental un- 
knowns consist of the seven generalized displace- 
ments and seven internal forces. The generalized 
displacements consist of the axial and transverse dis- 
placement components of the centroid u°, t>°, and 
w ° ; the three rotation components (jP x , 0°, and 0°; 
and the rate of twist of the beam 0 ° . The seven in- 
ternal forces consist of the axial force N x ; the two 
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transverse shear forces Q y and Q z \ the two bend- 
ing moments M y and M z ; the twisting moment M t ; 
and the bimoment . The sign convention for the 
generalized displacements and internal forces is given 
in figure 1. The fundamental equations of the thin- 
walled beam theory used herein are summarized in 
appendix A. In order to reduce the continuity re- 
quirements for all the generalized displacements to 
C° , a modified form of the Hellinger-Reissner mixed 
variational principle (Washizu 1974) is used in the de- 
velopment of the element characteristics. The mod- 
ification amounts to augmenting the functional of 
that principle by two terms: the Lagrange multiplier 
associated with the constraint condition relating <j ) ° 
and 0°, and a regularization term that is quadratic 
in the Lagrange multiplier. (See appendix A.) For 
a detailed discussion of the augmented (perturbed) 
Lagrangian formulation, see Fortin and Glowinski 
(1983). 

The shape functions used in approximating the 
internal forces and the Lagrange multiplier are se- 
lected to be the same but are different from those 
used in approximating the generalized displacements. 
Moreover, the continuity of both the internal forces 
and the Lagrange multiplier is not imposed at in- 
terelement boundaries. 

The finite element equations for each individual 
element can be cast in the following compact form: 



where {i/}, {A}, and {X} are the internal force 
parameters, Lagrange multiplier parameters, and 
nodal displacements, respectively; [F] is the 
matrix of linear flexibility coefficients; [ S ] is the 
strain-displacement matrix; [P] and [Q\ are matrices 
associated with the constraint condition and the reg- 
ularization term in the functional, respectively (see 
appendix A); [M] is the consistent mass matrix; £2 is 
the frequency of vibration; £ is a penalty parameter 
associated with the regularization term; and a dot 
(•) refers to a zero submatrix. The formulas for the 
elemental arrays [F], [5], [P], [Q], and [M] are given 
in appendix B. 

If the internal force parameters and the La- 
grange multiplier parameters are eliminated from 
equations (1), one obtains the following equation in 
the nodal displacements { X }: 


([S]'[F]- l [S]-e[0]tP] -1 W]‘){X} = n 2 [M]{*} (2) 


Comments on Perturbed Lagrangian 
Formulation 

The following three comments regarding the per- 
turbed Lagrangian formulation used herein seem to 
be in order: 

1. The governing finite element equations of the 
perturbed Lagrangian formulation (eqs. (1)) include 
both those of the Lagrange multiplier approach and 
the penalty method as special cases as follows: 

(1) By letting the penalty parameter e go to 
infinity, equations (1) reduce to those of the Lagrange 
multiplier approach. 

(2) By eliminating the Lagrange multiplier 
terms from equations (1), the resulting equations are 
identical to those of the penalty method. 

2. The perturbed Lagrangian formulation, in ad- 
dition to reducing the continuity requirements for 
the generalized displacements in the element devel- 
opment, alleviates two of the drawbacks of the La- 
grange multiplier approach and the penalty method, 
namely: 

(1) The regularization term in the functional 
results in replacing one of the zero diagonal blocks in 
the discrete equations of the Lagrange multiplier ap- 
proach by the diagonal matrix [P]/e in equations (1). 

(2) The constraint condition relating the ro- 
tation of the cross section <f> x and the rate of twist 
0° (see appendix A) is enforced exactly in the func- 
tional rather than approximately as in the penalty 
method. However, the presence of the regularization 
term results in replacing the constraint condition by 
the perturbed condition (on the element level): 

H{A} + [Q]'{*} = 0 (3) 

3. An important consideration in the perturbed 
Lagrangian formulation is the proper selection of the 
penalty parameter e. With the foregoing mixed mod- 
els, the penalty parameter can be chosen indepen- 
dently of the element size without adversely affect- 
ing the performance of the model. The accuracy of 
the solution increases with increasing value of the 
penalty parameter. However, for very large values 
of e, ill-conditioning of the equations occurs, thereby 
increasing the round-off errors. In the present study, 
for convenience, the penalty parameter e was chosen 
to be equal to the extensional stiffness of the beam 
EA. 

Numerical Studies 

To test and evaluate the performance of the fore- 
going mixed models, a large number of vibration 
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problems of thin-walled beams with open cross sec- 
tion have been solved by using these models. For 
each problem, comparison was made with solutions 
obtained by using two-dimensional plate/shell ele- 
ments in modeling the beam, and also with previ- 
ously published solutions whenever available. Herein, 
the results of two sets of free-vibration problems are 
discussed. The two problems are: (1) a straight can- 
tilever beam with symmetric and unsymmetric chan- 
nel cross sections; and (2) a semicircular beam with 
a Z-shaped section (Z-section) and clamped ends. A 
classical Rayleigh-Ritz solution along with test re- 
sults for the first problem are given by Gere and Lin 
(1958). Also, a finite element solution for the same 
problem based on displacement models is given by 
Mei (1970). 

In each case, the structure was analyzed us- 
ing (1) the mixed one-dimensional models developed 
herein, and (2) two-dimensional plate/shell models 
based on a Sanders-Budiansky type of shell theory 
with the effects of transverse shear deformation and 
rotary inertia included. (See Noor and Andersen 
1982 and Noor and Peters 1983 for a description 
of the finite element model used.) In the one- 
dimensional models, Lagrangian interpolation func- 
tions were used for approximating each of the gen- 
eralized displacements, internal forces, and Lagrange 
multiplier. The polynomial functions for the internal 
forces and the Lagrange multiplier were one degree 
lower than those of the generalized displacements. 
The characteristics of the finite element models in 
the numerical studies are summarized in table I. In 
all the solutions presented herein, only the displace- 
ment boundary conditions are enforced and the force 
boundary conditions are not applied. Typical re- 
sults are presented in figures 2 through 7 for the 
cantilever beams with channel sections and in fig- 
ures 8 through 12 for the semicircular beams with 
Z-sections. Also, table II lists the minimum six fre- 
quencies obtained by the thin-walled beam and two- 
dimensional models for the unsymmetric channel sec- 
tion and Z-section. 

Cantilever Beams With Channel Cross 
Sections 

The first problem set considered is that of the 
cantilever beams with the symmetric and the unsym- 
metric channel cross sections shown in figure 2. The 
vibration mode shapes for the beam with symmet- 
ric channel section, obtained by the two-dimensional 
model, are shown in figure 3. The ratios of the 
strain energies associated with normal stresses and 
torsional shear stresses to the total strain energy, 
U n /U and U ts /U, for the first six vibration modes 
are shown in figure 4. Note that U n is associated 


with N x , M y , M z , and and U i8 is associated 
with M t . Except for the second and sixth modes 
of the beam with symmetric channel section, all the 
other modes exhibit strong flexural-torsional cou- 
pling. This strong coupling cannot be predicted by 
ordinary beam models in which B w is neglected. For 
all the six modes, the transverse shear strain energy 
was found to be negligible (less than 0.6 percent of 
the total energy). 

Figures 5 and 6 give an indication of the accu- 
racy and convergence of the first six frequencies ob- 
tained by using the M3-2 and M4-3 models. (See ta- 
ble I.) The standard for comparison was taken to be 
the two-dimensional-model solution obtained using 
a 20 x 12 grid of M9-4 elements (20 elements along 
the beam length). As can be seen from figures 5 
and 6, the frequencies predicted by the mixed mod- 
els were highly accurate and converged rapidly with 
the increase in the number of elements. Thus, the 
converged frequencies obtained by the thin-walled 
beam model were close to those obtained by the two- 
dimensional model. The maximum error in the low- 
est six frequencies obtained by the converged thin- 
walled beam models was less than 1 percent. 

The effect of the magnitude of the penalty pa- 
rameter on the accuracy of the lowest six frequencies 
obtained by the mixed models is depicted in figure 7. 
As can be seen, the accuracy of the frequencies ob- 
tained by the mixed models is fairly insensitive to the 
choice of e in the range of e/EA from 10“ 3 to 10 2 . 
The upper limit of e is a function of the number of 
significant digits in the computer. 

Semicircular Beams With Z-Sections 

The second problem is that of the clamped semi- 
circular beams with the Z-sections of equal and un- 
equal flanges. (See fig. 8.) Henceforth, the equal 
and unequal flange sections will be referred to as 
the “symmetric and unsymmetric Z-sections.” Be- 
cause of the spanwise symmetry of the structure, only 
half the beam was modeled. The vibration mode 
shapes associated with the lowest six frequencies for 
the beam with symmetric Z-section are shown in fig- 
ure 9. These mode shapes are obtained by the two- 
dimensional model. The ratios of the strain energies 
U n /U and U ts /U obtained by the thin-walled beam 
model for the six vibration modes are shown in fig- 
ure 10. Again, for all six modes the transverse shear 
strain energy was found to be less than 0.6 percent 
of the total energy. 

An indication of the accuracy and convergence of 
the frequencies obtained by the M3-2 and M4-3 mod- 
els for the beams with symmetric and unsymmetric 
sections is given in figures 11 and 12. Although the 
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frequencies predicted by the thin-walled beam mod- 
els converge rapidly with increasing the number of 
elements, the errors in the converged frequencies ob- 
tained by the thin-walled beam model are greater 
than those for the straight beam described in the 
previous subsection. The maximum error in the low- 
est six frequencies for the beam with unsymmetric 
Z-section was of the order of 7.5 percent. (See fig. 12 
and table II.) The accuracy degradation for the beam 
with unsymmetric Z-section may be attributed to the 
cross-sectional distortions that are more pronounced 
in the mode shapes of the semicircular beam than in 
those of the straight beam. 

Concluding Remarks 

Simple, mixed finite element models are devel- 
oped for the free-vibration analysis of curved, thin- 
walled beams with arbitrary open cross section. The 
analytical formulation is based on a linear, Vlasov- 
type, thin-walled beam theory with the effects of 
flexural-torsional coupling, transverse shear deforma- 
tion, and rotary inertia included. The fundamental 
unknowns consist of seven internal forces and seven 
generalized displacements of the beam. The element 
characteristic arrays are obtained by using a modified 
form of the Hellinger-Reissner mixed variational prin- 
ciple. The modification consists of augmenting the 
functional of that principle by two terms: (1) the La- 
grange multiplier associated with the constraint con- 
dition relating the rotation of the cross section and 
the twist degrees of freedom, and (2) a regularization 


term that is quadratic in the Lagrange multiplier. 
Only C° continuity is required for the generalized 
displacements. The internal forces and the Lagrange 
multiplier are allowed to be discontinuous at interele- 
ment boundaries. 

The high accuracy and effectiveness of the mixed 
models developed are demonstrated by means of nu- 
merical examples of thin- walled beams with symmet- 
ric and unsymmetric cross sections. The standard of 
comparison in these examples is taken to be the solu- 
tions obtained by using two-dimensional plate/shell 
models for the beams. The accuracy of the mixed 
models developed is dependent on the cross-sectional 
distortions during the beam deformations. As the 
cross-sectional distortions increase, the degradation 
of accuracy becomes more pronounced. 
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Appendix A 

Fundamental Equations of Thin-Walled Beam Theory Used in Present Study 

The fundamental equations of the linear, Vlasov-type theory of curved thin-walled beams are given in this 
appendix. A right-handed orthogonal coordinate system is used with the x-axis passing through the centroids 
of the cross sections. (See fig. 1.) The beam is assumed to be curved in one direction only (in the £ 2 -piane). 

Displacement Assumptions 

Based on the assumption that the projection of each cross section on a plane normal to the initial centroidal 
axes does not distort during deformation, the displacement field in the plane of the cross section (yz- plane) is 
represented by 



where u°, v 0 , and up are the axial and transverse displacement components at y = z = 0; (p®, <t>y> and <f > ® 
are the rotation components about the coordinate axes; 0 0 is the rate of twist of the beam; and u is the 
sectorial coordinate (warping of the cross section for a unit rate of twist). The seven generalized displacement 
parameters u°, v°, w 4%, <j>% 4>2> an< i 0° are functions of x only. 

Strain Assumptions 


The following expressions are used for the three nonzero components of the strain field in the plane of the 
cross section: 


£x=e%- y*y + ZK°z - 
Ixy = 1% ~ ZK t 
Ixz = Ixz + V K t 


(A2) 


where e° is the extensional strain of the centerline; k® and /c® are the curvature changes in the y- and z- 
directions, respectively; /c® is the twist; and 7^y and 7 ^ are the transverse shear strains. The strain parameters 
e®, Ky, Kg, 7 j 2 , 7° y , , and ip° are functions of x only and can be expressed in terms of the displacement and 

rotation components as follows: 


£ 0 = 3u° + — 

x R 

Ky ~ 0(P Z ^ 

K 0 g = djPy 
Ixy = dv ° ~ <t>°z 


Ixz = + dw ° + <t>y 

V»° = d0° 


(A3) 


where d = d/dx and R is the radius of curvature of the centerline of the beam. Also, the following constraint 
condition is used to relate 6° and <f>®: 

d<t% - 0° = 0 (A4) 
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Constitutive Relations 

The relations between the internal forces and the strain components are given by 



where A is the cross-sectional area; is the first sectorial moment of the cross section; J y , I z , and I yz are the 
second moments of the cross section (moments and product of inertia); J is the Saint- Venant torsion constant; 
/<*,, Ivy i and Iuj Z are the second sectorial moments of the cross section (sectorial moments of inertia); E and 
G are the Young’s and shear moduli of the material, respectively; N x is the axial force; M y and M z are the 
bending moments; B w is the bimoment; Q y and Q z are the transverse shearing forces; and Mt is the twisting 
moment. The definition of the sectorial properties of the cross section is given by Friberg (1985), Heins and 
Wang (1978), Vlasov (1961), and Zbirohowski-Koscia (1967). 

Variational Functional 

The functional used in the element development is given by 

* = 7THR + J o A ( d<j>x ~ 0 °) dx ~^J Q (^) 2 dx 


where A is the Lagrange multiplier, e is a penalty parameter, l is the length of the element, and 7 Thr is the 
functional of the Hellinger-Reissner mixed variational principle. The expression for 7 Tjjr is 

*HR= f l {V-U c + K)dx (A7) 

where 0 



(A8) 


(A9) 
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(A10) 


K = f fi 2 [a [( U y + ( V y + (™ 0 ) 2 ] + (/, + /,) (<j>° x ) 2 + Iy (4$f + u (*!) 2 - 2Iy*m 

+ 2 (S„U° - I^yty + Iu,Z<t>i) 0° + U0 0 ) 2 } 

where p is the mass density of the material. Note in equations (A5) through (A10) that y and z are centroidal 
coordinates. (See fig. 1.) The sectorial properties of the cross section are evaluated using the centroid as the 
pole and the origin of the contour coordinates f as the lower free end of the cross section. (See fig. 1.) A 
FORTRAN program for evaluating the sectorial properties is listed by Coyette (1987). 
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Appendix B 

Formulas for Coefficients in Governing Equations for Individual Elements 

The explicit forms of the elemental arrays [F], [5], [P], [Q], and [ M ] are given in this appendix. For 
convenience, each of these arrays is partitioned into blocks according to contributions from individual nodes, 
internal forces, or Lagrange multiplier parameters. The expressions of the typical partitions (or blocks) are 
given in table BI. Note that the order of internal force parameters is N x , M z , M v , B w , Q y , Q z , and M t ; the 
order of the generalized displacement coefficients is u ° , v° , w°, 4>° x , <^J> $!> an< ^ • 

In table BI, N* and N J are the shape functions for both the internal forces and the Lagrange multiplier; 
N* represents the shape functions for the generalized displacements; R is the radius of curvature; m is the 
number of displacement nodes in the element; [/] is the identity submatrix; s is the number of parameters used 
in approximating each of the internal forces and the Lagrange multiplier; l is the length of the element; and 
d = d/dx. The range of the indices i and j is 1 to s, and the range of the indices i and j is 1 to m. The dots 
in the matrices refer to zero terms. 


Table BI. Explicit Form of Typical Partitions of the Arrays [F], [5], [P], [Q], and [M] 



\ 

I 

j 

I 
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Table I. Characteristics of Mixed Finite Element Models Used in Numerical Studies 


Model 

Number of 
displacement 
nodes 

Number of parameters 
per internal force 
(or Lagrange multiplier) 

Number of 
quadrature 
points 

Designation 

One-dimensional 

3 

2 

3 

M3-2 


4 

3 

4 

M4-3 

Two-dimensional 

9 

4 

9 

M9-4 


Table II. Minimum Frequencies Obtained by Two-Dimensional 
and Thin- Walled Beam (Converged) Models for the Cantilever 
Beam With Unsymmetric Channel Section (Fig. 2) and the 
Semicircular Beam With Unsymmetric Z-Section (Fig. 8) 


Mode 

Minimum frequencies 
for cantilever beam, Hz 

Minimum frequencies 
for semicircular beam, Hz 

^2d 

^converged 

n 2 d 

^converged 

i 

11.64 

11.66 

12.14 

12.48 

2 

19.64 

19.67 

23.12 

23.08 

3 

40.94 

41.09 

44.14 

42.71 

4 

58.52 

58.97 

50.90 

50.88 

5 

103.6 

103.9 

92.71 

91.42 

6 

113.9 

114.6 

100.2 

92.71 
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Figure 1. Characteristics of thin- walled beam and sign convention for internal forces and displacements. 




Figure 2. Cantilever beams with channel cross sections used in present study. The subscript s refers to 
shear center. 




frequencies obtained by two-dimensional model. 
















I 
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Figure 6. Accuracy and convergence of frequencies obtained by mixed models for cantilever beam with 
unsymmetric channel section shown in figure 2. 







E = 7.17 x 10 10 Pa 
\> = 0.3 

P = 2768 kg/m 3 
R = 1.0 m 
b = 5.08 x 10' 2 m 
t = 6.0 x 10‘^m 

Figure 8. Semicircular beams with Z-sections used in present study. The subscript s refers to shear center. 


Symmetric Unsymmetric 

section section 
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Q 1 = 10.70 (11.10) 
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Figure 9. Vibration mode shapes for semicircular beam with symmetric Z-section shown in figure 8. 
Numbers in parentheses are converged thin-walled beam frequencies (in hertz); other numbers are 
frequencies obtained by two-dimensional model. 
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